#(GPL, copyleft, https://www.gnu.org/licenses/gpl-3.0.de.html)
This code shows model comparision and confirmatory factor analyses.
Copyright (C) <2021>  <Walter, B.; Zehtner, R. I. & Hermann, A.>
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

#This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

#You should have received a copy of the GNU General Public License
along with this program.  If not, see 
https://www.gnu.org/licenses/ <https://www.gnu.org/licenses/>

#================================================================

#analyses were computed with R version 4.0.4 (R Core Team, 2019) 

#To run this program you need to specify the path of data as well as the data frames
#Note: Items 19 and 28 were deleted before running the CFAs

#================================================================

rm(list = ls())
#setwd('D:/FEQ_GR/Study1/Analyses')

library(tidyverse)
library(foreign)
library(psych)
library(psychTools)
library(lattice)

library(lavaan) # for CFA
library(semTools) # for additional functions in SEM
library(semPlot) # for path diagram

library(nonnest2) #for model comparison
library(MBESS) #for CI of reliability estimates

# library(gmodels)
# library(summarytools)
# library(RColorBrewer)
# library(haven)



# ++++++++++++++++++++++
# create data frames
# ++++++++++++++++++++++
df1 <- read.spss( "FEQ-GR_subsample1.sav", use.value.labels=FALSE, to.data.frame=TRUE )
names(df1)
feq_sample_1 <- df1[36:73]

df2 <- read.spss( "FEQ-GR_subsample2.sav", use.value.labels=FALSE, to.data.frame=TRUE )
names(df2)
feq_sample_2 <- df2[36:73]

feq_total <- bind_rows(feq_sample_1, feq_sample_2)

str(feq_total)

# ++++++++++++++++++++++
# Model comparison 
# ++++++++++++++++++++++

#### 4 factors: PD, PS, ND, NS

fa4_model_1 <-"
fa4_scalePD =~ FEQ_01 + FEQ_06 + FEQ_16 + FEQ_17 + FEQ_18 + FEQ_23 + FEQ_26 + FEQ_33 + 
               FEQ_39
               
fa4_scalePS =~ FEQ_02 + FEQ_03 + FEQ_13 + FEQ_21 + FEQ_22 + FEQ_30 + FEQ_31 + 
               FEQ_35 + FEQ_38 + FEQ_40             

fa4_scaleND =~ FEQ_04 + FEQ_05 + FEQ_07 + FEQ_09 + FEQ_11 + FEQ_12 + FEQ_24 + FEQ_27 + 
              FEQ_36 + FEQ_37
              
fa4_scaleNS =~ FEQ_08 + FEQ_10 + FEQ_14 + FEQ_15 + FEQ_20 + FEQ_25 + 
              FEQ_29 + FEQ_32 + FEQ_34
"


### 3 factors: P, NS, ND

fa3_model_1 <-"
fa3_scaleP =~ FEQ_01 + FEQ_06 + FEQ_16 + FEQ_17 + FEQ_18 + FEQ_23 + FEQ_26 + FEQ_33 + 
               FEQ_39 + FEQ_02 + FEQ_03 + FEQ_13 + FEQ_21 + FEQ_22 + FEQ_30 + FEQ_31 + 
               FEQ_35 + FEQ_38 + FEQ_40

fa3_scaleND =~ FEQ_04 + FEQ_05 + FEQ_07 + FEQ_09 + FEQ_11 + FEQ_12 + FEQ_24 + FEQ_27 + 
              FEQ_36 + FEQ_37
              
fa3_scaleNS =~ FEQ_08 + FEQ_10 + FEQ_14 + FEQ_15 + FEQ_20 + FEQ_25 + 
              FEQ_29 + FEQ_32 + FEQ_34
"


### 2 factors: P, N

fa2_model_1 <-"
fa2_scaleP =~ FEQ_01 + FEQ_06 + FEQ_16 + FEQ_17 + FEQ_18 + FEQ_23 + FEQ_26 + FEQ_33 + 
               FEQ_39 + FEQ_02 + FEQ_03 + FEQ_13 + FEQ_21 + FEQ_22 + FEQ_30 + FEQ_31 + 
               FEQ_35 + FEQ_38 + FEQ_40

fa2_scaleN =~ FEQ_04 + FEQ_05 + FEQ_07 + FEQ_09 + FEQ_11 + FEQ_12 + FEQ_24 + FEQ_27 + 
              FEQ_36 + FEQ_37 + FEQ_08 + FEQ_10 + FEQ_14 + FEQ_15 + FEQ_20 + FEQ_25 + 
              FEQ_29 + FEQ_32 + FEQ_34
"



#fit2 = 2 factor model, fit3 = 3 factor model, fit4 = 4 factor model

fit2 <- cfa(fa2_model_1, data = feq_sample_1, estimator = "ML", mimic = "Mplus")
fit3 <- cfa(fa3_model_1, data = feq_sample_1, estimator = "ML", mimic = "Mplus")
fit4 <- cfa(fa4_model_1, data = feq_sample_1, estimator = "ML", mimic = "Mplus")


###fit4: covariance matrix of latent variables is not positive definite

summary(fit4, fit.measures=T, standardized=T)

###### comparision only between fit2 and fit3 possible

vuongtest(fit2, fit3) 


# ++++++++++++++++++++++
####### CFA subsample 1
# 2 factors
# ++++++++++++++++++++++

fa2_model_1 <-"
fa2_scaleP =~ FEQ_01 + FEQ_06 + FEQ_16 + FEQ_17 + FEQ_18 + FEQ_23 + FEQ_26 + FEQ_33 + 
               FEQ_39 + FEQ_02 + FEQ_03 + FEQ_13 + FEQ_21 + FEQ_22 + FEQ_30 + FEQ_31 + 
               FEQ_35 + FEQ_38 + FEQ_40

fa2_scaleN =~ FEQ_04 + FEQ_05 + FEQ_07 + FEQ_09 + FEQ_11 + FEQ_12 + FEQ_24 + FEQ_27 + 
              FEQ_36 + FEQ_37 + FEQ_08 + FEQ_10 + FEQ_14 + FEQ_15 + FEQ_20 + FEQ_25 + 
              FEQ_29 + FEQ_32 + FEQ_34
"

s1.fa2.cfa1 = cfa(fa2_model_1, data = feq_sample_1, estimator = "MLR", mimic = "Mplus")

summary(s1.fa2.cfa1, fit.measures = T, standardized = T)

##reliability

s1.fa2.cfa1.rel = reliability(s1.fa2.cfa1)
print(s1.fa2.cfa1.rel, digits = 3)

##CI reliability
feq_sample_1_fa2_scaleP_items <- feq_sample_1 %>% select(FEQ_01, FEQ_06,
                                                         FEQ_16, FEQ_17, FEQ_18, FEQ_23, FEQ_26, FEQ_33,
                                                         FEQ_39,
                                                         FEQ_02, FEQ_03, FEQ_13, FEQ_21, FEQ_22, FEQ_30, FEQ_31,
                                                         FEQ_35,
                                                         FEQ_38, FEQ_40)

feq_sample_1_fa2_scaleN_items <- feq_sample_1 %>% select(FEQ_04, FEQ_05,
                                                         FEQ_07, FEQ_09, FEQ_11, FEQ_12, FEQ_24, FEQ_27,
                                                         FEQ_36,
                                                         FEQ_37, FEQ_08, FEQ_10, FEQ_14, FEQ_15, FEQ_20, FEQ_25,
                                                         FEQ_29,
                                                         FEQ_32, FEQ_34)


set.seed(42)
ci.reliability(data = feq_sample_1_fa2_scaleP_items,
               type="hierarchical", conf.level = 0.95, interval.type="bca", B=1000)

set.seed(42)
ci.reliability(data = feq_sample_1_fa2_scaleN_items,
               type="hierarchical", conf.level = 0.95, interval.type="bca", B=1000)


# ++++++++++++++++++++++
#######CFA subsample 1
# 3 Faktoren: P, NS, ND
# ++++++++++++++++++++++

fa3_model_1 <-"
fa3_scaleP =~ FEQ_01 + FEQ_06 + FEQ_16 + FEQ_17 + FEQ_18 + FEQ_23 + FEQ_26 + FEQ_33 + 
               FEQ_39 + FEQ_02 + FEQ_03 + FEQ_13 + FEQ_21 + FEQ_22 + FEQ_30 + FEQ_31 + 
               FEQ_35 + FEQ_38 + FEQ_40

fa3_scaleND =~ FEQ_04 + FEQ_05 + FEQ_07 + FEQ_09 + FEQ_11 + FEQ_12 + FEQ_24 + FEQ_27 + 
              FEQ_36 + FEQ_37
              
fa3_scaleNS =~ FEQ_08 + FEQ_10 + FEQ_14 + FEQ_15 + FEQ_20 + FEQ_25 + 
              FEQ_29 + FEQ_32 + FEQ_34
"

s1.fa3.cfa1 = cfa(fa3_model_1, data = feq_sample_1, estimator = "MLR", mimic = "Mplus")

summary(s1.fa3.cfa1, fit.measures = T, standardized = T)

##reliability
s1.fa3.cfa1.rel = reliability(s1.fa3.cfa1)
print(s1.fa3.cfa1.rel, digits = 3)

##CI reliability
feq_sample_1_fa3_scaleP_items <- feq_sample_1 %>% select(FEQ_01, FEQ_06, 
                                                         FEQ_16, FEQ_17, FEQ_18, FEQ_23, FEQ_26, FEQ_33,
                                                         FEQ_39, 
                                                         FEQ_02, FEQ_03, FEQ_13, FEQ_21, FEQ_22, FEQ_30, FEQ_31,
                                                         FEQ_35, 
                                                         FEQ_38, FEQ_40)

feq_sample_1_fa3_scaleND_items <- feq_sample_1 %>% select(FEQ_04, FEQ_05, FEQ_07, FEQ_09, FEQ_11, FEQ_12, FEQ_24, FEQ_27,
                                                          FEQ_36, FEQ_37)

feq_sample_1_fa3_scaleNS_items <- feq_sample_1 %>% select(FEQ_08, FEQ_10, FEQ_14, FEQ_15, FEQ_20, FEQ_25, 
                                                          FEQ_29, FEQ_32, FEQ_34)


set.seed(42)
ci.reliability(data = feq_sample_1_fa3_scaleP_items, 
               type="hierarchical", conf.level = 0.95, interval.type="bca", B=1000)

set.seed(42)
ci.reliability(data = feq_sample_1_fa3_scaleND_items, 
               type="hierarchical", conf.level = 0.95, interval.type="bca", B=1000)

set.seed(42)
ci.reliability(data = feq_sample_1_fa3_scaleNS_items, 
               type="hierarchical", conf.level = 0.95, interval.type="bca", B=1000)



# ++++++++++++++++++++++
#######CFA subsample 2
#3 Scales: 
##P: -item 13, 
##ND: no deletion of items, 
##NS: -item 29
# ++++++++++++++++++++++

s1.items_f3 <- select(feq_sample_2, -FEQ_13, -FEQ_29)

fa3_model_2 <-"
fa3_scaleP =~ FEQ_01 + FEQ_06 + FEQ_16 + FEQ_17 + FEQ_18 + FEQ_23 + FEQ_26 + FEQ_33 + 
  FEQ_39 + FEQ_02 + FEQ_03 + FEQ_21 + FEQ_22 + FEQ_30 + FEQ_31 + 
  FEQ_35 + FEQ_38 + FEQ_40

fa3_scaleND =~ FEQ_04 + FEQ_05 + FEQ_07 + FEQ_09 + FEQ_11 + FEQ_12 + FEQ_24 + FEQ_27 + 
  FEQ_36 + FEQ_37

fa3_scaleNS =~ FEQ_08 + FEQ_10 + FEQ_14 + FEQ_15 + FEQ_20 + FEQ_25 + 
  FEQ_32 + FEQ_34
"

s1.fa3.cfa2 = cfa(fa3_model_2, data = s1.items_f3, estimator = "MLR", mimic = "Mplus")

summary(s1.fa3.cfa2, fit.measures = T, standardized = T)

##reliability
s1.fa3.cfa2.rel = reliability(s1.fa3.cfa2)
print(s1.fa3.cfa2.rel, digits = 3)

##CI reliability
feq_sample_2_fa3_scaleP_items <- feq_sample_2 %>% select(FEQ_01, FEQ_06, 
                                                         FEQ_16, FEQ_17, FEQ_18, FEQ_23, FEQ_26, FEQ_33,
                                                         FEQ_39, 
                                                         FEQ_02, FEQ_03, FEQ_21, FEQ_22, FEQ_30, FEQ_31,
                                                         FEQ_35, 
                                                         FEQ_38, FEQ_40)

feq_sample_2_fa3_scaleND_items <- feq_sample_2 %>% select(FEQ_04, FEQ_05, FEQ_07, FEQ_09, FEQ_11, FEQ_12, FEQ_24, FEQ_27,
                                                          FEQ_36, FEQ_37)
feq_sample_2_fa3_scaleNS_items <- feq_sample_2 %>% select(FEQ_08, FEQ_10, FEQ_14, FEQ_15, FEQ_20, FEQ_25, FEQ_32, FEQ_34)


set.seed(42)
ci.reliability(data = feq_sample_2_fa3_scaleP_items, 
               type="hierarchical", conf.level = 0.95, interval.type="bca", B=1000)

set.seed(42)
ci.reliability(data = feq_sample_2_fa3_scaleND_items, 
               type="hierarchical", conf.level = 0.95, interval.type="bca", B=1000)

set.seed(42)
ci.reliability(data = feq_sample_1_fa3_scaleNS_items, 
               type="hierarchical", conf.level = 0.95, interval.type="bca", B=1000)




